Librairies et jeu de données

library(rmarkdown)
library(dplyr)
library(tidyverse) # pour  %>%
library(kableExtra) # kbl
library(knitr) # kable
library(ggplot2)
library(plotly) # graphes dynamiques
library(MASS)
library(corrplot) # corrplot
library(GGally) # ggpairs
library(gridExtra) # co-plot pour ggplot
library(ggfortify) # autoplot
library(car) # Tests type II
library(rsample)   # create training/test sets 
library(corrplot)

1. Motivations

1.1.Corrélation/causalité

Corrélation : La corrélation mesure la relation statistique entre deux variables, indiquant si et dans quelle mesure elles varient ensemble.

Causalité : La causalité désigne la relation de cause à effet où un événement (cause) entraîne un autre événement (effet), établissant une relation de dépendance directe.

⚠️ La corrélation ne prouve pas la causalité.


1.2. Niveau d’education vs Chomage

Ici on observe une corrélation entre deux variables “Plus d’éducation = moins de chômage”, mais cela ne signifie pas nécessairement qu’il y a une relation de cause à effet (causalité) entre ces variables :

  • Corrélation Spurieuse : Des facteurs non mesurés tels que l’économie, la politique ou le social peuvent influencer à la fois l’éducation et le chômage sans qu’il y ait un lien direct entre eux.

  • Direction de la Causalité Inverse : Il est possible que ce soit une économie prospère qui encourage davantage l’éducation plutôt que l’éducation seule qui réduise le chômage.

  • Facteurs Confondants : Des variables telles que la situation économique globale ou les politiques gouvernementales peuvent influencer à la fois l’éducation et le chômage, créant ainsi une corrélation apparente.

  • Causalité Complexes : La relation entre l’éducation et le chômage peut être complexe, impliquant des interactions multiples et non linéaires avec d’autres facteurs, ce qui rend difficile l’attribution d’une causalité basée uniquement sur la corrélation observée.

En résumé, bien qu’il puisse exister une corrélation entre le niveau d’éducation et le taux de chômage, d’autres facteurs doivent être pris en compte pour déterminer s’il existe réellement une relation de cause à effet entre ces variables. Des études plus approfondies, des analyses de causalité et la prise en compte de variables confondantes sont nécessaires pour établir des liens de causalité robustes.

1.3. Nombre de prix Nobel vs Consommation de chocolat

Dans ce cas on observe une corrélation linéaire entre deux variables telles que le nombre de prix Nobel dans un pays et la consommation de chocolat, il est peu probable que cela implique une relation de causalité directe entre ces variables. Cette situation illustre le principe selon lequel la corrélation ne prouve pas la causalité.

  • Corrélation Fortuite : La corrélation entre la consommation de chocolat et les prix Nobel peut être due au hasard.
  • Corrélation Induite : Une forte consommation de chocolat et un nombre élevé de prix Nobel peuvent être liés par d’autres facteurs comme le niveau de développement économique.
  • Absence de Mécanisme Causal : Il n’y a pas de lien logique direct entre la consommation de chocolat et les réalisations des lauréats de prix Nobel.

En résumé, bien que la corrélation puisse être observée entre le nombre de prix Nobel dans un pays et la consommation de chocolat, il est important de ne pas conclure à une relation de causalité directe sur la seule base de la corrélation. Des études plus approfondies, une analyse des mécanismes causaux potentiels et la prise en compte d’autres variables sont nécessaires pour évaluer la nature de la relation entre ces variables.

1.4. L’exemple de l’Ozone

Ce jeu de données (Air Breizh, Summer 2001) contient 112 observations et 13 variables

ozone1 <- read.table("ozone1.txt",header=TRUE, sep="", dec=",")
paged_table(ozone1)

Corrigeons la déclaration erronnée des variables:

#library(dplyr)
ozone1$Date <- as.Date(as.factor(ozone1$Date),format = "%Y%m%d")
ozone1$vent<-as.factor(ozone1$vent)
ozone1$pluie<-as.factor(ozone1$pluie)
ozone1 <- ozone1 %>% mutate(across(where(is.character), as.numeric))
ozone1 <- ozone1 %>% mutate(across(where(is.integer), as.numeric))
paged_table(ozone1)

💡 Bon à savoir !

  • Pour changer la nature en date, il faut d’abord transformer en facteur, puis préciser le format.

  • Ici, across est utilisé avec where(is.character) pour cibler uniquement les colonnes qui sont de type caractère (chr) et appliquer as.numeric à ces colonnes.


But Prédire la concentration maximum d’ozone dans l’air (maxO3). Interessons nous dans un premier temps uniquement à la Températures relevées à 12h (T12). Trouver une fonction \(f\) telle que

maxO3\(\approx f\)(T12)

avec

  • \(X=\)T12 Températures relevées à 12h.
  • \(Y=\)maxO3 : Concentration maximum d’ozone dans l’air.

Comment choisir \(f\)

  • Fonctions affines (laquelle?)
  • Fonctions polynômes (degré?)
  • Spline (degré?)
  • Noyaux (lequel?)

\(\Rightarrow\) Trouver la fonction qui minimise la somme des carrés des distances de la fonction aux points (\(T12_i,maxO3_i\)).

2. Premiers Modèles de régression linéaire

2.1. Une seule variable explicative

Intéressons-nous au modèle linéaire, c’est-à-dire que

maxO3\(\approx\) fonction affine de T12 \(\Leftrightarrow\) \(Y=f(X)+\epsilon=\beta_0+\beta_1X+\epsilon\)

  • \(X=\)T12 : variable réponse, variable à expliquer,..
  • \(Y=\)maxO3 : variable explicative, régresseur, prédicteur, covariable,..
  • \(\epsilon\) : erreurs, résidu théorique
  • \(\beta_0\) : ordonnée à l’origine ou intercept
  • \(\beta_1\): pente ou poids de la variable \(X\)

Problème \((\beta_O,\beta_1)\) sont inconnus!

Solution les estimer à partir de nos données, de notre échantillon de \(n\) observations \((X_i,Y_i)_{i=1,\cdots,n}\).

Modèle linéaire simple:

\[Y_i=f(X_i)+\epsilon_i=\beta_0+\beta_1X_i+\epsilon_i,\qquad i=1,\cdots,n\]

[Hypothèses:]

  • [P1] Les \(\epsilon_i\) sont centrées
  • [P2] Les \(\epsilon_i\) ne sont pas corrélées
  • [P4] Les \(\epsilon_i\) ont une variance (niveau de bruit) constante \(\sigma^2>0\) (homoscédasticité)
  • [P4] Les \(\epsilon_i\) sont gaussiennes.

On peut résumer ces hypothèses par \(\epsilon_i\overset{iid}{\sim}\mathscr N(0,\sigma^2)\)

[Prédiction le modèle estimé:] \[\widehat Y_i=\widehat f\big(X_i\big)= \widehat\beta_0+\widehat\beta_1X_{i}\]


Pour trouver \((\widehat\beta_0,\widehat\beta_1)\), on minimise la somme des erreurs aux carrés à partir du modèle \[Y_i=f(X_i)+\epsilon_i=\beta_0+\beta_1X_i+\epsilon_i,\qquad i=1,\cdots,n\] C’est à dire \((\widehat\beta_0,\widehat\beta_1)\) sont les arguments qui minimise la somme des carrés résiduels \[ (\widehat\beta_0,\widehat\beta_1)=\underset{(\beta_0,\beta_1)}{\arg\min}\sum_{i=1}^n(Y_i- f(X_i))^2=\underset{(\beta_0,\beta_1)}{\arg\min}\sum_{i=1}^n(Y_i-\beta_0-\beta_1X_i)^2=\underset{(\beta_0,\beta_1)}{\arg\min}\sum_{i=1}^n\epsilon_i^2 \] On appelle cet estimateur, l’estimateur des moindres carrés, qui est également dans ce cadre l’estimateur du maximum de vraisemblance. Avec R:

mod<-lm(maxO3~T12,data=ozone1)
mod$coefficients%>%kbl(col.names=c("Estimation"))
Estimation
(Intercept) -27.419636
T12 5.468685


Traçons cette droite, appellée la droite des moindres carrées.

ggplot(ozone1)+ aes(x = T12, y =maxO3)+ geom_point(col='red', size = 0.5)+ geom_smooth(method = "lm",se = FALSE)

2.2. Deux variables explicatives (numérique+facteur)

Si maintenant on considère nos observations en prenant en compte le facteur ventet que l’on souhaite faire régression différente pour chaque modalité.

p1<-ggplot(ozone1)+ aes(x = T12, y =maxO3)+geom_point(size = 0.5)
p2<-p1+  geom_smooth(method = "lm",se = FALSE)
p3<-ggplot(ozone1)+ aes(x = T12, y =maxO3, color=vent)+geom_point(size = 0.5)
p4<- p3+guides(color = FALSE)
p5<-p3+  geom_smooth(method = "lm",se = FALSE)
subplot(p1,p2,p4,p5,nrows = 2)

Alors il faut définir autant de modèle linéaire du type

\[Y_{i}=f(X_i)+\epsilon_i=\beta_{0}+\beta_{1}X_{i}+\epsilon_{i}, i=1,\cdots,n\] qu’il y a de modalités \(k\).

Modèle linéaire (Forme régulière):

\[Y_{ik}=f(X_{ik})+\epsilon_{ik}=\beta_{0k}+\beta_{1k}X_{ik}+\epsilon_{ik},\qquad i=1,\cdots,n_k, \ k=1,\dots,4\] où

  • \(\beta_{0k}\) est l’intercept pour les observations admettant la modalité numéro \(k\) pour le facteur vent .
  • \(\beta_{1k}\) est la pente (=poids de la variable \(X\)) pour les observations admettant la modalité numéro \(k\) pour le facteur vent .


Mais on pourrait imaginer qu’il existe un intercept et une pente de reférence identiques à tous

  • \(\beta_{0k}=\mu+\alpha_k\): où \(\mu\) est l’intercept de référence et les \(\big(\alpha_k\big)_k\) l’impact du facteur vent.
  • \(\beta_{1k}=b+c_k\): où \(b\) est la pente de référence et les \((c_k)_k\) l’interaction entre la variable \(X\)=T12et le facteur vent.

Modèle linéaire (Forme singulière):

Pour tout \(i=1,\cdots,n_k,\) et tout \(k=1,\dots,4\) \[Y_{ik}=f(X_{ik})+\epsilon_{ik}=\mu+\alpha_k+(b+c_k)X_{ik}+\epsilon_{ik}\]

⚠️ Ce modèle n’est pas identifiable,il existe une infinité de solutions. Pour le rendre identifiable, il faut ajouter des contraintes.


Par défaut R ajoute la contrainte \(\alpha_1=c_1=0\) pour “exhiber” une solution.

lm(maxO3~T12+vent+T12:vent, data=ozone1)

ou bien

lm(maxO3~T12*vent, data=ozone1)

mod_defaut<-lm(maxO3~T12*vent, data=ozone1)
mod_defaut$coefficients%>%kbl(col.names=c("Estimation"))
Estimation
(Intercept) 8.735714
T12 4.092281
ventNord -24.902135
ventOuest -57.366714
ventSud -43.195199
T12:ventNord 1.006046
T12:ventOuest 2.232719
T12:ventSud 1.680646


  • \(\alpha_1=\)ventEst n’apparait pas car R a pris par défault \(\widehat\alpha_1=0\)
  • \(b_1=\)T12:ventEst n’apparait pas car R a pris par défault \(\widehat b_1=0\)

Retrouver le modèle régulier il suffit sous R d’indiquer que \(\mu=b=0\)

lm(maxO3~-1+T12+vent, data=ozone1)

mod_reg<-lm(maxO3~-1+T12+vent, data=ozone1)
library(kableExtra)



mod_reg$coefficients%>%kbl(col.names=c("Estimation"))
Estimation
T12 5.479838
ventEst -24.107753
ventNord -23.821255
ventOuest -30.814974
ventSud -27.504905


  • Ìntercept n’apparait pas car R a pris par défault \(\widehat\mu=0\)
  • \(b=\)T12 n’apparait pas car R a pris par défault \(\widehat b=0\)

💡 Bon à savoir !

Quel que soit le modèle et les contraintes choisis, la prédiction sera la même ; seule l’interprétation des paramètres peut différer.


3. Modèle de régression linéaire multivarié (ML)

3.1. Déclaration du modèle

But Prédire la concentration maximum d’ozone dans l’air (maxO3) en fonction des variables de notre jeu de données à partir d’un modèle linéaire.

  • Eviction de la date On exclue tout d’abord la variable Date
ozone1bis<-ozone1[,-c(1)]
  • train/test On sépare notre jeu de donnée en 2 (80%-20%) pour garder des données pour évaluer la qualité de notre modèle.
set.seed(2024)
ozone_split    <- rsample::initial_split(ozone1bis, prop = 0.8)
ozone<- rsample::training(ozone_split)
ozone_test <- rsample::testing(ozone_split)
  • \(Y=\)maxO3 : Concentration maximum d’ozone dans l’air.
  • \(X^{(1)}=\) T9 Températures relevées à 9h.
  • \(X^{(2)}=\) T12 Températures relevées à 12h.
  • \(X^{(3)}=\) T15 Températures relevées à 15h.
  • \(X^{(4=)}=\) Ne9 Variable synthétique nebuleuse à 9h.
  • \(X^{(5)}=\) Ne12 Variable synthétique nebuleuse à 12h.
  • \(X^{(6)}=\) Ne15 Variable synthétique nebuleuse à 15h.
  • \(X^{(7)}=\) Vx9 Variable synthétique vent à 9h.
  • \(X^{(8)}=\) Vx12 Variable synthétique vent à 12h.
  • \(X^{(9)}=\) Vx15 Variable synthétique ventà 15h.
  • \(X^{(10)}=\) maxO3v .
  • \(X^{(11)}=\) vent Facteur vent à 4 modalités (Est, Nord, Ouest, Sud).
  • \(X^{(12)}=\) pluie Facteur pluie à 2 modalités (Pluie, Sec).

maxO3\(=f\)(T12,T9,T12,T15,Ne9,Ne12,Ne15,Vx9,Vx12,Vx15,maxO3v,vent,pluie)\(+\epsilon\)

\[ Y=f\big(X^{(1)},\cdots,X^{(12)}\big)+\epsilon \]

où les mêmes hypothèses sont faites:

[Hypothèses:]

  • [P1] Les erreurs \(\epsilon\) sont centrées
  • [P2] Les erreurs \(\epsilon\) ne sont pas corrélées
  • [P4] Les erreurs \(\epsilon\) ont une variance (niveau de bruit) constante \(\sigma^2>0\) (homoscédasticité)
  • [P4] Les erreurs \(\epsilon\)sont gaussiennes.

Prédiction \(\widehat Y_i=\widehat f\big(X^{(1)},\cdots,X^{(12)}\big)\)


Si on ne considère aucune interactions alors sous R le modèle par défaut est le modèle singulier et s’écrit

lm(maxO3~T9+T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v+vent+pluie,data=ozone)

ou bien plus rapidement

lm(maxO3~.,data=ozone)

Déclarons notre modèle complet mod_completet affichons la sortie summary(mod_complet)

mod_complet<-lm(maxO3~., data=ozone)
summary(mod_complet)
## 
## Call:
## lm(formula = maxO3 ~ ., data = ozone)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.518  -8.506  -0.731   7.450  44.092 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  22.2487    18.1291   1.227   0.2236    
## T9            1.5073     1.4334   1.052   0.2964    
## T12           0.6172     1.8575   0.332   0.7406    
## T15           0.2322     1.4355   0.162   0.8720    
## Ne9          -2.6971     1.1829  -2.280   0.0255 *  
## Ne12          0.4995     1.7572   0.284   0.7770    
## Ne15         -0.6239     1.2956  -0.482   0.6316    
## Vx9           0.4761     1.0939   0.435   0.6647    
## Vx12          0.6873     1.4162   0.485   0.6289    
## Vx15          0.4899     1.0771   0.455   0.6505    
## maxO3v        0.3822     0.0775   4.931 4.87e-06 ***
## ventNord     -0.8189     7.8379  -0.104   0.9171    
## ventOuest     3.1013     9.4555   0.328   0.7438    
## ventSud       4.6715     8.4367   0.554   0.5814    
## pluieSec      2.9572     4.1185   0.718   0.4750    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.88 on 74 degrees of freedom
## Multiple R-squared:  0.7505, Adjusted R-squared:  0.7033 
## F-statistic:  15.9 on 14 and 74 DF,  p-value: < 2.2e-16

3.2. Interprétation de la sortie

Regardons de plus près


On observe 4 colonnes

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.2487    18.1291   1.227   0.2236    
T9            1.5073     1.4334   1.052   0.2964    
T12           0.6172     1.8575   0.332   0.7406    
T15           0.2322     1.4355   0.162   0.8720    
Ne9          -2.6971     1.1829  -2.280   0.0255 *  
Ne12          0.4995     1.7572   0.284   0.7770    
Ne15         -0.6239     1.2956  -0.482   0.6316    
Vx9           0.4761     1.0939   0.435   0.6647    
Vx12          0.6873     1.4162   0.485   0.6289    
Vx15          0.4899     1.0771   0.455   0.6505    
maxO3v        0.3822     0.0775   4.931 4.87e-06 ***
ventNord     -0.8189     7.8379  -0.104   0.9171    
ventOuest     3.1013     9.4555   0.328   0.7438    
ventSud       4.6715     8.4367   0.554   0.5814    
pluieSec      2.9572     4.1185   0.718   0.4750    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  • Estimate : correspond aux \(\widehat\beta_j\) l’estimation des paramètres \(\beta_j\) par maximum de vraissemblance/moindres carrés
  • Std. Error : Racine carré de l’erreur d’estimation des \(\beta_j\)
  • t value: Valeur de la statistique de test du test : \(H_0\, : \ \beta_j=0 \ \ vs \ H_1\, : \ \ \beta_j\neq 0\)
  • Pr(>|t|) : \(p\)-value du test (de Student) précédent.

⚠️ Attention !

Ne pas rejeter \(H_0\) \(\large{\neq}\) accepter \(H_1\). Ne signifie pas que la variable n’est pas pertinente (écarter une variable est ici prématuré).


💡 Bon à savoir ! Pour afficher uniquement ces 4 colonnes

mod_complet%>% summary()%>%coefficients()%>%kbl()


Dans la sortie du summary(mod_complet) nous avons également la ligne suivante :

F-statistic:  15.9 on 14 and 74 DF,  p-value: < 2.2e-16

Test de Fisher global: Teste la pertinence des variables dans leur ensemble (sous \(H_0\) modèle réduit à l’intercept)}.

  • F-statistic: 15.9 : Valeur de la statistique de test: \(H_0\, : \ \forall j: \ \beta_j=0 \ \ vs \ H_1\, : \overline{H_0}\)
  • p-value: < 2.2e-16 : \(p\)-value du test. Ici, on rejette \(H_0\), il existe au moins une varibale pertinente.

💡 Bon à savoir !

Le modèle réduit à l’intercept s’écrit

mod_0<-lm(maxO3~1,data=ozone)


Dans la sortie du summary(mod_complet) nous avons également les lignes suivantes :

Residual standard error: 14.88 on 74 degrees of freedom
Multiple R-squared:  0.7505,    Adjusted R-squared:  0.7033
  • Residual standard error : Estimation \(\widehat\sigma\) du niveau de bruit \(\sigma\) (écart type) défini dans les hypothèses du modèle. Ici \(\widehat\sigma=14.88\), très élevé
  • Multiple R-squared : coefficient determination \(R^2\)
  • Adjusted R-squared : coefficient determination ajusté \(R^2_a\)

Interprétation des coeff de determination: Ils sont compris entre 0 et 1. La valeur 0 indique que le modèle ne colle pas aux données, la valeur 1indique que le modèle parfaitement aux données.

⚠️ Attention ! Plus le nombre de variables explicatives est grand, plus le \(R^2\) est grand! Le coefficient \(R^2_a\) est moins impacté par ce nombre (=la dimension)


Pour finir, dans la sortie du summary(mod_complet) nous avons également les lignes suivantes :

Residuals:
    Min      1Q  Median      3Q     Max 
-48.518  -8.506  -0.731   7.450  44.092
  • Residuals : Les résidus estimés (les théoriques \(\epsilon_i\) sont inconnus) définis de la façon suivantes: \[\widehat\epsilon_i=Y_i-\widehat Y_i=Y_i-\widehat f\big(X^{(1)},\cdots,X^{(12)}\big)\]

📌 Pour valider le modèle (c’est-à-dire les hypothèses), il faudra faire une analyse des résidus. Mais avant de valider le modèle, avons-nous sélectionner le “meilleur” modèle ? Toutes les variables sont-elles pertinentes ?


3.3. Etude des corrélations

library(dplyr)
library(corrplot)
ozone_numerique<- ozone  %>%  select_if(is.numeric)

correlation_matrix<-cor(ozone_numerique)
corrplot(correlation_matrix, order = 'hclust',addrect = 3)

Les variables numériques sont très corrélées entre elles!

ggpairs(ozone_numerique)

plot_ly(y = ozone$maxO3, color = ozone$pluie,type = 'box') %>%
  layout(yaxis = list(title = 'max03'))
plot_ly(y = ozone$maxO3, color = ozone$vent,type = 'box') %>%
  layout(yaxis = list(title = 'max03'))

4. Sélection de modèles

4.1. Methodes pas à pas

Méthodes pas-à-pas:

On choisit un critère (\(R^2_a\), AIC (par défaut), BIC, \(C_p\)-Mallows,…) ET une des 3 méthodes pas-à-pas suivantes

  • Méthode Backward : On part du modèle de taille \(p\) variables explicatives et on enlève 1 à 1 les variables en comparant les modèles 2 à 2 selon un critère.

  • Méthode Forward : On part du modèle de taille \(1\) réduit à l’intercept et on ajoute 1 à 1 les variables en comparant les modèles 2 à 2 selon un critère.

  • Méthode Both: mixte de 2 méthodes. On part de l’intercept et on ajoute/enlève les variables 1 à 1 et on compare selon un critère.

Arrêt lorsque le critère n’est plus améliorable .


mod_0=lm(maxO3~1,data=ozone)

mod_forw<-step(mod_0, maxO3~T9+T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v+vent+pluie,data=ozone,trace=F,direction=c('forward'))
mod_forw
## 
## Call:
## lm(formula = maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
## 
## Coefficients:
## (Intercept)          T12       maxO3v          Ne9         Vx12  
##     20.6041       2.2690       0.3923      -2.5949       1.0888
mod_both<-step(mod_0,maxO3~T9+T12+T15+Ne9+Ne12+Ne15+Vx9+Vx12+Vx15+maxO3v+vent+pluie,data=ozone, trace=F,direction=c("both"))
mod_both
## 
## Call:
## lm(formula = maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
## 
## Coefficients:
## (Intercept)          T12       maxO3v          Ne9         Vx12  
##     20.6041       2.2690       0.3923      -2.5949       1.0888
mod_back<-step(mod_complet,~.,trace=F,data=ozone,direction=c("backward"))
mod_back
## 
## Call:
## lm(formula = maxO3 ~ T9 + Ne9 + Vx9 + maxO3v, data = ozone)
## 
## Coefficients:
## (Intercept)           T9          Ne9          Vx9       maxO3v  
##     24.4931       2.6017      -2.9336       1.6811       0.3845

Les 2 méthodes donne le même modèle, c’est un hasard. Mais comparer au modèle complet et celui restreint à l’intercept, quel est le meilleur model?

AIC(mod_complet,mod_0,mod_forw,mod_both,mod_back)

Notre modèle réduit à 4 variables + l’intercept est meilleur selon le critère AIC.

summary(mod_both)
## 
## Call:
## lm(formula = maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -54.128  -9.989  -0.907   8.630  43.605 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 20.60408   12.63631   1.631  0.10673    
## T12          2.26903    0.54497   4.164 7.55e-05 ***
## maxO3v       0.39230    0.06874   5.707 1.68e-07 ***
## Ne9         -2.59485    0.87285  -2.973  0.00385 ** 
## Vx12         1.08876    0.70288   1.549  0.12514    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.33 on 84 degrees of freedom
## Multiple R-squared:  0.7375, Adjusted R-squared:  0.725 
## F-statistic: 59.01 on 4 and 84 DF,  p-value: < 2.2e-16

Notons que si le \(R^2\) a diminué le \(R^2_a\) a augmenté. Les erreurs d’estimation ont également diminuées. Le modèle est meilleur.

4.2. Tests de Type I et II

Test de type I

anova(mod_complet)

Ligne à ligne un test est fait avec sous \(H_0\) un modèle avec \(d+1\) variables vs sous \(H_1\) un modèle avec \(d+2\). On ajoute les variables 1 à une, pour tester leur pertinence

  • Ligne 1: \(H_0\) intercept vs \(H_1\) intercept+T9
  • Ligne 2: \(H_0\) intercept+T9 vs \(H_1\) intercept+T9+T12
  • Ligne 3: \(H_0\) intercept+T9+T12 vs \(H_1\) intercept+T9+T12+T15
  • Ligne 4: \(H_0\) intercept+T9+T12+T15 vs \(H_1\) intercept+T9+T12+T15+Ne9
  • Ligne 5: \(H_0\) intercept+T9+T12+T15+Ne9 vs \(H_1\) intercept+T9+T12+T15+Ne9+Ne12
  • …

👀 Ici on garde les variables : T9, T12, Ne9, Vx9et maxO3v. Définissons ce modele

model_typeI<-lm(maxO3~T9+T12+Ne9+Vx9+maxO3,data=ozone)

Test de type II Ligne à ligne un test est fait avec sous \(H_0\) un modèle avec \(p\) variables vs sous \(H_1\) un modèle avec \(p+1\). On retire une variable pour tester sa pertinence

  • Ligne 1: \(H_0\) modele complet vs \(H_1\) modele complet sans T9
  • Ligne 2: \(H_0\) modele complet vs \(H_1\) modele complet sans T12
  • Ligne 3: \(H_0\) modele complet vs \(H_1\) modele complet sans T15
  • Ligne 4: \(H_0\) modele complet vs \(H_1\) modele complet sans Ne9
  • Ligne 5: \(H_0\) modele complet vs \(H_1\) modele complet sans Ne12
  • …
library(car)
Anova(mod_complet)

👀 Ici on garde les variables : Ne9 et maxO3v. Définissons ce modele

model_typeII<-lm(maxO3~Ne9+maxO3,data=ozone)

Comparons nos modeles

AIC(mod_complet,mod_0,mod_both,model_typeI,model_typeII)

Le modèle mod_both reste le meilleur suivant le critère AIC.

5. Validation du modèle

library(ggfortify)
autoplot(mod_both,1:3)

Hypothèse [P1]: Erreurs sont centrées/le modèle est linéaire.

  • Si dans le plot Residuals vs Fitted les résidus globalement uniformément répartis des deux côtés de 0. Si la ligne bleue est approximativement horizontale et à zéro (échelle!) alors [P1] est validée.

📌 Les fitted sont les réponses prédites \(\widehat Y_i\) par le modèle.

span class=“solution”>👀 Ici [P1] validée

Hypothèse [P2]: Variance est homoscédastique.

  • Si dans le plot Scale-Location, la courbe bleue est approximativement horizontale et à 1 (échelle!) et que les points uniformément répartis autour de celle-ci (Absence de tendance). alors [P2] est validée.

span class=“solution”>👀 Ici [P2] validée

  • Si il est difficile de valider alors faire test de Breush-Pagan pour décider. \(H_0\): “Homoscédasticité”
ncvTest(mod_both)$p
## [1] 0.006507032

span class=“solution”>👀 Ici \(p\)-value<5 % donc [P2] non validée

Hypothèse [P3]: Erreurs sont non corrélées.

  • Faire un test de Durbin-Watson pour décider. \(H_0\): “Non-corrélation”
set.seed(2024)
durbinWatsonTest(mod_both)$p
## [1] 0.264

span class=“solution”>👀 Ici \(p\)-value>5 % donc [P3] validée

Hypothèse [P4]: Erreurs sont gaussiennes.

  • Si dans le plot Normal Q-Q:n, les points sont alignés autour de la première bissectrice alors [P4] est validée.

  • Si il est difficile de valider alors faire test de Shapiro pour décider. \(H_0\): “Gaussien”

shapiro.test(residuals(mod_both))$p
## [1] 0.01259202

span class=“solution”>👀 Ici \(p\)-value>5 % donc [P4] non validée

6. Identification et gestion des outliers

6.1. Identification

Outliers à fort résidus:

Le modèle peut échoué à prédire correctement \(Y_i\) \(\Rightarrow\) grande différence entre \(|Y_i-\widehat Y_i|=|\epsilon_i|>2\). On peut identifier ces observations sur le Studentized-plot.


influenceIndexPlot(mod_both,vars="Studentized")

2 observations sont des outilers à fort résidus : 31et 37

Outliers à fort effet levier sur eux-même:

Si une observation avec grande influence sur sa propre estimation (\(\widehat Y_i\)), on dit qu’elle à un fort effet levier sur sa propre estimation : \(\Rightarrow\) la hat-value>0.5. On peut identifier ces observations sur le Hat-plot


influenceIndexPlot(mod_both,vars="hat")

Aucune observation n’a de fort impact sur sa propre estimation.

Outliers à fort effet levier sur le modèle:


Si une observation a une grande influence sur le modèle lui-même (\(\widehat \beta_j\)), on dit qu’elle à un fort effet levier sur le modèle : \(\Rightarrow\) si sa distance de Cook est grande. On peut identifier ces observations sur le Cook-plot

influenceIndexPlot(mod_both,vars="cook")

Aucune observations ne semble avoir une distance très grande. L’observation 31 a la distance de cook la plus grande.

Mais doit-on pour autant retirer ces points ? Tous les ouliers sont-ils à retirer ?

6.2. Gestion

Le critère de Bonferroni est une méthode d’ajustement pour tenir compte du fait que plusieurs tests sont effectués simultanément. Il contrôle le taux d’erreur de famille, c’est-à-dire la probabilité de commettre au moins une erreur de type I parmi tous les tests effectués.

Dans le contexte de influenceIndexPlot avec vars = "Bonf", les \(p\)-values associées à chaque observation sont ajustées selon le critère de Bonferroni. Les observations qui ont des \(p\)-values ajustées inférieures à un seuil spécifié sont considérées comme ayant une influence significative sur le modèle.

Lorsque vous utilisez influenceIndexPlot avec l’argument vars = "Bonf", cela signifie que vous souhaitez inclure des lignes dans le graphique des indices d’influence qui correspondent aux observations ayant une influence significative sur le modèle, selon le critère de Bonferroni.

Les outliers problématiques ont une \(p\)-value de Bonferroni \(<\) seuil On peut afficher le Bonferroni-plot

influenceIndexPlot(mod_both,vars="Bonf")

Il est important de noter que l’interprétation des p-values ajustées dépend du seuil alpha choisi. Un seuil alpha plus strict (par exemple, 0.01) conduira à moins d’observations étiquetées comme ayant une influence significative, tandis qu’un seuil alpha plus permissif (par exemple, 0.05) conduira à plus d’observations étiquetées.

En résumé, lorsque vous utilisez vars = "Bonf" avec influenceIndexPlot, les observations étiquetées ont des p-values ajustées selon le critère de Bonferroni qui sont inférieures au seuil alpha spécifié, ce qui indique une influence significative sur le modèle de régression. Adjuster le seuil alpha permet de contrôler la sensibilité du critère d’inclusion des observations dans le graphique des indices d’influence.

Mais le plus simple est d’acceder directement aux \(p\)-value des points problématiques. Pour afficher les \(p\)-\(values\) des points jugés douteux par R, il suffit d’écrire

outlierTest(mod_both)
##     rstudent unadjusted p-value Bonferroni p
## 31 -4.418044         2.9857e-05    0.0026573

Dans notre exemple: l’observation 31 a la \(p\)-value au sens de Bonferroni.

ozone[31,]

7. Prédiction et évaluation de notre modèle

mod_complet<-lm(maxO3~., data=ozone)
mod_F_compl<-lm(maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone)
mod_F_ss_outlier<-lm(maxO3 ~ T12 + maxO3v + Ne9 + Vx12, data = ozone[-c(31),])

7.1. Prediction

Ozone_actual<-ozone_test$maxO3
Results<-as.data.frame(Ozone_actual)
Results$Modele_complet<-predict(mod_complet,newdata=ozone_test)
Results$Modele_Final<-predict(mod_F_compl,newdata=ozone_test)
Results$Modele_Final_ss_outlier<-predict(mod_F_ss_outlier,newdata=ozone_test)

Results
library(ggplot2)

# Plot pour le modèle complet
p1 <- ggplot(Results, aes(x = Ozone_actual, y = Modele_complet)) +
  geom_point(color = "#31AA41") +geom_smooth(method = "lm", se = FALSE,color = "#31AA41") +ylab("Modele complet") 

# Plot pour le modèle avec outlier
p2 <- ggplot(Results, aes(x = Ozone_actual, y = Modele_Final)) +geom_point(color = "#AA0012") +geom_smooth(color = "#AA0012",method = "lm", se = FALSE) + ylab("Modele avec outlier") 

# Plot pour le modèle sans outlier
p3 <- ggplot(Results, aes(x = Ozone_actual, y = Modele_Final_ss_outlier))+
  geom_point(color = "#4200AA") +geom_smooth(color = "#4200AA",method = "lm", se = FALSE) +ylab("Modele ss outlier") 

# Affichage des plots côte à côte
grid.arrange(p1, p2, p3, nrow = 2)

7.2. Evaluation

Calculons l’erreur moyenne faite par chacun de ces modèles sur ces données non vues

library(Metrics)
Erreur_mod_complet<-rmse(ozone$maxO3,predict(mod_complet,newdata=ozone))
Erreur<-as.data.frame(Erreur_mod_complet)
Erreur$Modele_Final<-rmse(ozone$maxO3,predict(mod_F_compl,newdata=ozone))
Erreur$Modele_Final_ss_outlier<-rmse(ozone$maxO3,predict(mod_F_ss_outlier,newdata=ozone))

Erreur_mod_complet<-rmse(ozone_test$maxO3,predict(mod_complet,newdata=ozone_test))
Erreur_test<-as.data.frame(Erreur_mod_complet)
Erreur_test$Modele_Final<-rmse(ozone_test$maxO3,predict(mod_F_compl,newdata=ozone_test))
Erreur_test$Modele_Final_ss_outlier<-rmse(ozone_test$maxO3,predict(mod_F_ss_outlier,newdata=ozone_test))

rbind.data.frame(train=Erreur,test=Erreur_test)

Pour aller plus loin

Random forest

library(randomForest)
library(Metrics)
# Créer un modèle de Random Forest
mod_RF <- randomForest(maxO3~.,data =ozone)
#c(rmse(ozone$maxO3),predict(mod_RF,newdata=ozone))
Erreur_mod_RF<-rmse(ozone$maxO3,predict(mod_RF,newdata=ozone))
Erreur_mod_RF_test<-rmse(ozone_test$maxO3,predict(mod_RF,newdata=ozone_test))
paste0('RF train:   ',Erreur_mod_RF,'   RF test:  ' ,Erreur_mod_RF_test)
## [1] "RF train:   6.45839847004582   RF test:  14.6830507786875"

👀 ICI les RF font de l’overfitting, le rmsesur le train et le test sont trop différent.

📌 Utilisation des RF pour faire de la selection. Utile lorsque le nombre de variables potentiellement explicatives est important.


library(vip)
vip(mod_RF,num_features = 14, bar = FALSE)

Si on prend le même nombre de variables que dansq le ML, c’est-à-dire 4, les RF n’améliorent pas notre prédiction

# Obtenir l'importance des caractéristiques
B<-vip(mod_RF,num_features = 4)
B<-B$data$Variable
ozone_RF<- ozone[,c('maxO3',B )]
mod_LM_RF<-lm(maxO3~.,data = ozone_RF)

Erreur_mod_LM_RF<-rmse(ozone$maxO3,predict(mod_LM_RF,newdata=ozone))
Erreur_mod_LM_RF_test<-rmse(ozone_test$maxO3,predict(mod_RF,newdata=ozone_test))
paste0('LM_RF train:   ',Erreur_mod_LM_RF,'   LM_RF test:  ' ,Erreur_mod_LM_RF_test)
## [1] "LM_RF train:   15.6072400107394   LM_RF test:  14.6830507786875"

Régression polynomiale

Reprenons pour la suite le jeu de données cars de la librairie ‘MASS’ qui contient

50 Données et 2 variables : speed:la vitesse en miles par heure et dist la distance nécessaire à l’arrêt d’un véhicule.

library(rmarkdown)
library(MASS)
data_cars <- as.data.frame(cars)
paged_table(data_cars)
library(ggplot2)
library(gridExtra)

plot_p<-ggplot(cars, aes(x = speed, y = dist)) + geom_point( size=0.85)
plot1 <- plot_p +  geom_smooth( color="#FF9900", method = "lm" )+  labs(title = "Régression linéaire")

plot2 <-plot_p+ geom_smooth(color="#4169E1", method = "lm",  formula = y ~ poly(x, 2, raw = TRUE)) +labs(title = "Polynôme de degré 2")

grid.arrange(plot1, plot2, nrow = 1)

Définissons ces modèles

mod<- lm(dist~speed,data=cars)
mod2<- lm(dist~poly(speed, 2),data=cars)

Regardons le plot des résidus

library(ggfortify)
autoplot(mod, 1:4)

autoplot(mod2, 1:4)

La courbe bleue à une forme parabolique. Essayons alors d’exprimer disten fonction d’un polynôme d’ordre 2

summary(mod)
## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12
summary(mod2)
## 
## Call:
## lm(formula = dist ~ poly(speed, 2), data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.720  -9.184  -3.188   4.628  45.152 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       42.980      2.146  20.026  < 2e-16 ***
## poly(speed, 2)1  145.552     15.176   9.591 1.21e-12 ***
## poly(speed, 2)2   22.996     15.176   1.515    0.136    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.18 on 47 degrees of freedom
## Multiple R-squared:  0.6673, Adjusted R-squared:  0.6532 
## F-statistic: 47.14 on 2 and 47 DF,  p-value: 5.852e-12

Regardons le plot des résidus pour valider le modèle

AIC(mod,mod2)
modele_lineaire<-rmse(cars$dist,mod$fitted.values)
Erreur<-as.data.frame(modele_lineaire)
Erreur$modele_polynomial<-rmse(cars$dist,mod2$fitted.values)
rownames(Erreur) <- "Erreur"
Erreur